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Abstract 

This paper introduces a new metric and mean on the set of positive 
semidefinite matrices of fixed-rank. The proposed metric is derived 
from a well-chosen Riemannian quotient geometry that generalizes 
the reductive geometry of the positive cone and the associated nat- 
ural metric. The resulting Riemannian space has strong geometrical 
properties: it is geodesically complete, and the metric is invariant 
with respect to all transformations that preserve angles (orthogonal 
transformations, scalings, and pseudoinversion) . A meaningful ap- 
proximation of the associated Riemannian distance is proposed, that 
can be efficiently numerically computed via a simple algorithm based 
on SVD. The induced mean preserves the rank, possesses the most 
desirable characteristics of a geometric mean, and is easy to compute. 



1 Introduction 

Positive definite matrices have become fundamental computational objects 
in many areas of engineering and applied mathematics. They appear as 
covariance matrices in statistics, as elements of the search space in convex and 
semidefinite programming, as kernels in machine learning, and as diffusion 
tensors in medical imaging, to cite a few. Computing with positive definite 
matrices involves approximations, interpolation, filtering, and estimation, 
leading to a realm of metric-based algorithms. In the recent years, it has been 
increasingly recognized that the Euclidean distance does not suit best the 
set of positive definite matrices — the positive symmetric cone Pn — and that 
working with the proper geometry does matter in computational problems. 
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The (non- Euclidean) natural metric of the positive cone [13] proceeds 
from the rich quotient geometry of this set and its structure of reductive 
homogeneous space (see Appendix 7.2 for more details). The resulting "nat- 
ural" metric is invariant to the action (by congruence) of the general linear 
group, a feature at the core of many desirable properties. The associated Rie- 
mannian mean can be called geometric rather than arithmetic. Recent con- 
tributions that have advocated the use of this metric in applications include 
[26, 14, 34] for tensor computing in medical imaging and [7] for radar process- 
ing. More theoretical results can be found in [21, 4, 6, 28]. Most notably, the 
natural metric corresponds to the Fisher information metric for the multivari- 
ate normal distribution [10, 31]. The comprehensive paper [32] uses this result 
to derive an intrinsic Cramer-Rao bound for the Gaussian covariance matrix 
estimation problem. Finally, the natural metric coincides with the metric 
defined by the natural self-concordant logarithmically homogeneous barrier 
(— log det A) on the symmetric cone (which is a convex set) in optimization 
[9, 22]. In particular, the interest in short-step methods relies on the property 
that Pn with its natural metric is geodesically complete, i.e., every maximal 
geodesic is defined for all i e M, and thus the boundary can not be reached 
in finite time. 

Because matrix algorithms tend to be applied to computational problems 
of ever-increasing size, they need to be adapted to remain tractable. Typical 
matrix computations (hke SVD, EVD, QR factorization, etc.) require 0{n^) 
operations for a positive definite matrix of size n, which limits their use in 
large-scale problems. A sensible remedy is to work with low-rank approxi- 
mations instead. A rank p approximation of a positive definite matrix can be 
factored as ^4 = ZZ^, where the matrix Z e R"^^ is of much reduced size if 
p <^ n, leading to a reduction of the numerical cost of typical matrix opera- 
tions from O(n^) to 0{np'^). If p is kept to a moderate value, the complexity 
of the resulting algorithms grows only linearly with the size of the problem. 

The natural metric developed for positive definite matrices is only valid 
for full-rank matrices. The goal of this paper is to extend the natural metric 
of the positive cone to the set S+(p, n) of symmetric positive semidefinite 
matrices of fixed-rank p < n. The set S"'"(p, n) admits a quotient geome- 
try that generalizes the quotient geometry of the cone in a way that pre- 
serves most of the desirable properties of the cone geometry. Motivated 
by the natural metric of the positive cone, the proposed geometry differs 
from the quotient geometry recently proposed [3]. The resulting "natural" 
metric — which, to the best of the authors' knowledge, has not appeared in 
the literature previously — preserves not all but a remarkably large number 
of invariance properties of the natural metric in the cone. More precisely, 
it is invariant under all transformations that preserve angles, that is, rota- 
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tions, scalings, and pseudo-inversion. Endowed with this metric S"'"(p, n) is 
geodesically complete. Moreover, we propose a meaningful and numerically 
efficient approximation of the Riemannian distance. The induced notion of 
mean is is shown to present all the desirable properties of a geometric mean. 

The proposed natural metric on S"'"(p,n) is viewed as an important step 
to generalize several existing algorithms for positive definite matrices to the 
semidefinite case. This applies not only to the computational problems men- 
tioned previously, but also to a growing realm of matrix nearness problems 
based on the use of Bregman divergences. We mention in particular the recent 
paper [11] that leaves as an open question the characterization of different 
types of projections to compute distances onto important sets of matrices, 
such as the positive semidefinite cone. 

The rest of the paper is organized as follows. In section 2 we summarize 
the existing work on the positive definite cone Pn, concentrating on the fea- 
tures most relevant for the paper. To gain insight on the main issues faced 
when extending the geometry of Pn to S"'"(p, n), section 3 focuses on the sim- 
plest case of 2 X 2 matrices of rank 1 . A metric extending the natural metric 
on the cone is derived using polar decomposition. We provide geometrical 
and physical justifications for the associated invariance, distance, and mean 
properties. 

Building upon the polar decomposition of vectors in the plane, we develop 
in section 4 a quotient geometry for S'^(p, n). The resulting natural metric 
decomposes as the sum of the natural metric on the cone and the standard 
metric of the Grassman manifold. We prove that S"'"(p, n) endowed with this 
metric is a Riemannian manifold. 

Section 5 is devoted to the Riemannian distance associated to the metric. 
Because the explicit calculation of geodesies is out of reach, we construct 
special curves that approximate the geodesies. The length of those curves 
can be calculated by a SVD based algorithm of numerical complexity O(np^). 
It provides a meaningful notion of closeness in S"'"(p,n) which inherits the 
invariance properties of the Riemannian distance. The singularities of this 
measure are also characterized. As an aside, the manifold S'^(p, n) is proved 
to be geodesically complete. 

Section 6 provides a new definition of geometric mean between two ma- 
trices of S"'"(p, n) based on the measure of closeness. We prove that the mean 
preserves the rank, we argue it deserves the appellation "geometric", and 
we prove it is a generalization to rank-deficient positive symmetric matrices 
of the mean associated with the natural metric on the positive cone. From 
a computational viewpoint, a main advantage of the metric is that, after a 
suitable SVD, the computations of the closeness and mean decouple into two 
separate problems: computing a Riemannian distance and mean 1. in the 
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(lower dimension) cone Pp and 2. in the Grassman manifold of subspaces of 
dimension p in M". 

Concluding remarks and perspectives are discussed in section 7. 



• Pn is the set of symmetric positive definite n x n matrices. 

• S+(p, n) is the set of symmetric positive semidefinite nxn matrices of 
rank p <n. We will only use this notation in the case p < n. 

• Gl(n) is the general linear group, that is, the set of invertible nxn 
matrices. 

• M"^^ is the set of full rank n x p matrices. 

• Vn,p = 0(n)/0(n — p) is the Stiefel manifold; i.e., the set of n x p 
matrices with orthonormal columns: U^U — I p. 

• Gr(p, n) is the Grassman manifold, that is, the set of p-dimensional 
subspaces of W^. It can be represented by the equivalence classes 



• Sym(n) is the vector space of symmetric nxn matrices. 

• diag(Ai, . . . , A„) is the nxn matrix with the Aj's on its diagonal. / = 
diag(l, . . . , 1) is the identity matrix. 

• range (^4) is the subspace of spanned by the columns of ^4 e M'^^'^. 

• TxM. is the tangent space to the manifold M. at X. 

2 Riemannian distances and geometric means 
on the symmetric cone 

The geometry of the n-dimensional symmetric cone Pn has been well-studied 
in the literature. This section reviews some of its relevant features in view 
of the main developments of the present paper. Given a matrix A e Pn, a 
starting point is the matrix factorization 



1.1 



Notation 



Vn,p/0(p). 



A = ZZ'^ = {UR){URf = UR^U^ = i?'^ 



(1) 
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where A e P^,, Z e Gl(n), i?, i?' e P„, [/ e 0(n). The right and left polar 

decompositions Z = UR = R'U are unique [15], but the factorization A = 
ZZ^ is unaffected by orthogonal transformations Z ^ ZD with O G 0(n). 
The matrix equalities (1) underline the quotient geometry of the cone Pn- 

P„ ^ Gl(n)/0(n) ^ (0(n) x Pn)/0(n). (2) 

The characterization (2) encodes the rich geometry of Pn as a reductive 
homogeneous space, as discussed in [17, 32]. The most relevant consequence 
of this feature to the present context is the existence of a Gl(n)-invariant 
metric on the manifold Gl(n)/0(n). This metric is called the natural metric 
on the symmetric cone [13]. Up to a scaling factor it is also known as the 
affine-invariant metric [26], the Siegel metric in symplectic geometry, and it 
coincides with the metric given by the Fisher information matrix for Gaussian 
covariance matrix estimation [32]. We briefly summarize how it is derived 
and its main properties in the present context. 

In view of (1), Gl(n) has a transitive action on Pn via congruence 

A ^ LAL^, (3) 

where L e Gl(n) and any matrix A is brought back to the identity matrix 
choosing L = A^^l"^ with A^l"^ defined by the polar factor K in (1). Like- 
wise, any tangent vector X G Sym(n) at identity / G Pn can be transported 
to a tangent vector A}I'^XA!'I'^ G T^P„. At identity /, the Gl(n)-invariant 
metric is defined as the usual scalar product g^^{Xx^X-^ — Tr (XiXl") = 
Tr (X1X2). The invariance of the metrics then implies g^^{A}!'^XxA^!'^ , Al'I'^X^A^^'^) = 
Tr (X1X2), which can only be satisfied with the definition 

^^"(L*!, A) = Tr {p^A-^D^A-^) (4) 

at any arbitrary A G Pn- The invariance of the metric (4) has direct implica- 
tions on the expression of the geodesies and the accompanying Riemannian 
distance. The exponential map at / is the usual matrix exponential 

00 

expf"X = expX = ^(1/A;!)X^ 

A;=0 

The Frobenius norm is the geodesic length (i(expX, /) = hence 

the formula at identity 

d(A/) = ||iogA||^. 

Invariance of the metric again extends the characterization of geodesies at 
arbitrary A G Pn [21, 32]: 

exp^"(tX) = exp(a-^/'XA-i/2)Ai/2, t > 0, 
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and the corresponding geodesic distance 



dp„{A,B) = d{A-'/'BA-'/',I) = \\log{A-'/'BA-'/')y = /J] W(A,), 



where are the generahzed eigenvalues of the pencil A — XB, i.e., the roots 
of det{AB~^ — XI). Note that the distance is invariant with respect to ma- 
trix inversion (^4, B) i— > {A~^, B~^) because log (A/j) is invariant to inversion 
Afc ^ A^\ 

The geodesic characterization provides a closed-form expression of the 
Riemannian (Karcher) mean of two matrices A,B^ The geodesic A(t) 
linking A and B is 

A{t) = exp^-{tX) = A^/^exp{t\og{A-^/^BA-^/^))A^/^, 

where A-^/^x^-i/s = log(A-^/^BA-'^/^) e Sym(n). The midpoint is ob- 
tained for t = 1/2: 

AoB^ ^V2(^-l/25^-l/2^1/2^1/2_ (g) 

For a comprehensive treatment of the geometric versus arithmetic means 
of positive definite matrices, the reader is referred to [4, 21]. Prom a geomet- 
ric viewpoint, there arc numerous reasons to prefer the natural metric (and 
associated distance and mean) on the symmetric cone to the fiat metric asso- 
ciated to the distance ||/1 — i^Hp. Indeed the symmetric cone is not a vector 
space. The flat metric can still be used because the set is convex. However, 
many applications treating the space of covariance matrices as a vector space 
can yield to degraded algorithm performances [32] . A further undesirable fea- 
ture of the flat metric is that is not a geodesically complete space (see 
section 5) since the geodesic A + t{B — A) is not a positive matrix for all t. In 
contrast, it becomes gcodesicaUy complete with the natural metric. A practi- 
cal consequence is that the natural metric is well-suited to short-step interior 
point methods in Pn. Remarkably, it coincides with the metric deflned by 
the natural self- concordant logarithmically homogeneous barrier (—log det 
^4) on the symmetric cone (which is a convex set) in optimization [9, 22]. 
Purthermore, the invariance to the group action (3) (implying in particular 
invariance with respect to inversion) is desirable in numerous applications 
(see, e.g., [6, 26, 32, 7]). In particular, if A is a covariance matrix E(^;s^), 
the action corresponds to a change of basis z i— Lz. As a consequence the 
natural metric is well-suited to intrinsic estimation algorithms for covariance 
matrix estimation; see [32]. In [26, 14, 6], the authors argue that invariance 
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with respect to inversion is adapted to the physics of diffusion tensors related 
to medical imaging. Finally, the associated (geometric) mean has numerous 
desirable properties developed in section 6. It is useful to keep in mind that 
the geometric mean coincides with the classical geometric mean \/ab when 
a,b e Pi. 

3 Extending the metric: A geometric insight 
in the plane 

In order to inherit some of the nice properties of the natural distance on 
when the matrices are not full-rank, we seek a distance which is invariant to 
the Gl(n) group action. For tutorial purposes wc start by considering 2x2 
matrices of rank 1: S+(l,2). For any L e Gl(2) and A,B e S+(l,2) one 
wants 

d{A,B) = d{LAL^,LBL^). 
This objective is, however, too ambitious. Taking 

o)' ""^C 0' """^ "^^(o °)' 

A is unchanged by the transformation whereas LBL^ = ( ^ ^2 ) ■ Thus if the 
distance is continuous in the matrix elements (for example it is associated to 
a metric), we have in the limit d{A^ B) — d{A, >!) = as e — > 0. This proves 
that d is not a distance since A^ B. 

Mimicking the matrix decomposition (1), write A G S^(l,2): 

A — XX — ur u , 

where u — (cos(^), sin(^)) e C is a unit vector and {r,9) is the po- 
lar representation of x. Without loss of generahty 9 is equated to ^ -|- jtt, 
j G Z, since x and —x correspond to the same A. Thus S"'"(l,2) can be 
equated to the space x MP-*^, where WF^ is the real projective space of 
dimension 1 (lines of M^). The group action (3) corresponds to a change of 
basis for the vector x ^ Lx. We have proved in the last paragraph that 
there is no distance which is invariant under every change of basis for x. 
Nevertheless one can define a distance which is invariant by scaling and or- 
thogonal transformation, i.e., G = {/iR : {fi,R) G M x 0(2)}. A sensible 
metric is ds"^ = dO"^ + k {dr/ry = dO"^ + k c?(logr)^ with A; > 0. Let us find 
the distance associated to this metric. Let A = xixj and B = X2xi . Let 
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{ri,9i) denote their polar coordinates for i — 1,2. It is always possible to 
have \02 — Oi\ < 7r/2 possibly replacing 9i by ± tt. The distance is then 

4+(i^2)(A B) = 102 -ei\^ + k \\og{n/r2)\\ k>0. (7) 

where the first term penalizes the distance between the subspaces range (A) = 
span(a;i) and range(i?)=span(a;2). It ensures the invariance by rotation/orthogonal 
transformation (which only affects 9) while the second term ensures the in- 
variance by scaling (which only affects r). Note that when range(A)=range(i?) 
the induced distance on the 1-dimensional subspace corresponds to the nat- 
ural distance (5) between matrices of Pi. The metric above rewrites in the 
general form 

ds"^ = Tr {du^du) + k d{\og rf, k> 0. (8) 

A symmetry-based justification The choice of the metric (8) can be 
derived from necessary conditions imposed by the desired symmetries (in- 

variances). Let A = xixj , B = X2X2 with xi, X2 G M^. We seek an invariant 
distance between A and B. This is equivalent to seeking an invariant distance 
between xi and X2] i.e., a scalar G-invariant function on M^xM^, satisfying the 
conditions required to be a distance. If H is any subgroup of G, the distance 
must be in particular if-invariant. Let H — x SO (2) C G. The identifi- 
cation M-l = H can be obtained via polar coordinates (r, 6) G x = H. 
A standard result (see, e.g., [23, 8]) is that every ff-invariant scalar func- 
tion of {xi,X2) & H X H {in particular any distance) is a function of the 
scalar invariants {r2/ri,92 — 9i); i.e., the coordinates of Xi^ * X2, where 
* is the group multiplication of H. Thus any invariant distance d writes 
d{xi,X2) = h{r-2/ri,92 — ^1), and since d is a. distance it must be symmetric 
d{xi,X2) = d{x2,Xi). Thus necessarily d{xi,X2) = /(|log(r2/ri)|, ~ ^i|) 
with / positive, monotone in each of its arguments, and only /(O, 0) is equal 
to 0. It proves that the distance (7) is a prototype for every G-invariant 
distance d{A, B). 

A physical justification The invariance properties of the distance (7) 
are meaningful from a physical viewpoint. Suppose Z = ]K{zz'^) is the co- 
variance matrix of a stochastic variable ^ G (for instance the position 
of an object) with zero mean such that every realization of 2; is on a line. 
Then Z G S'^(l,2). The distance (7) between the two covariances Z\,Z2 
of two independent variables Z\,Z2 (dispersing in two different directions) is 
well-defined. It does not depend on any nontrivial way on the choice of mea- 
surement units, e.g., feet versus meters, as well as on the orientation of the 
frame chosen, e.g., the first axis is pointing north or south. 
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Moreover, suppose the measurements are noisy. For instance Zi — diag(4, e^) 
and Z2 = diag(e^, 1), where the term e ^ 1 is the amphtude of the covariance 
of the noise. These two matrices belong to P2, and the geometric mean (6) 
is diag(2e, e). The smallness of the noise ruins the mean, which no longer 
reflects the physical interpolation between the two processes. Indeed in the 
degenerate case e ^ 0, the mean becomes the null matrix. In contrast, the 
midpoint in the sense of the distance (7) between the rank 1 approximations 
of the two matrices (i.e., e = 0) is (11) with eigenvalues (0,2), which is an 
image of the covariance matrix of the midpoint. The distance between the 
matrices and their mean is independent of the choice of units and orientation, 
and is hardly affected by noise. 

4 A new Riemannian metric on the symmet- 
ric semidefinite cone 

Let S+(p, n) be the set of positive semidefinite matrices of fixed-rank p < n. 
One can prove (analogously to the last section) in the general case that is 
impossible to find a distance between matrices of S"'"(p, n) which would be 
invariant to the transformation (3) for p arbitrary small. Nevertheless one 
can construct a metric generalizing (8) and a distance generalizing (7). The 
Grassmanian distance is a multidimensional generalization of the angular 
distance d9, and the natural metric on the symmetric cone (4) is a multidi- 
mensional generalization of dlogr = r~^dr. The resulting metric is invariant 
to orthogonal transformations and scalings. Locally, when the ranges of two 
infinitely close matrices coincide, the induced metric on the corresponding 
subspace reduces to the natural metric on the symmetric cone Pp. In partic- 
ular, when p = n the metric coincides with the natural metric on Pn. 

Mimicking the developments of section 2, we start from the matrix equal- 
ities 

A = ZZ'^ = {UR){URf = UR^U'^, 

where ^ e S+(p,n), Z e WI'^p, ReV^.Ue Vn,p. Right multiphcation of 
Z by an orthogonal matrix Z 1— > ZO, O e 0(p), does not affect the product 
A — ZZ^ . Consider the following 0(p) group action: 

R ^ O'^RO e Pp, 

U^UO e Vn,p. 

The representation A — UR^U^ with {U, R^) G Vn,p x Pp is thus univocal up 
to the equivalence relation (C/, R^) = {UO, O^R^O) for any O e 0(p). Thus 
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the set S''"(p, n) admits a quotient manifold representation 



S+(p,n)-(Vn,pXPp)/0(p). 



(A dimension checking yields that the dimension of S+(p, n) is dim(Vn,p x Pp)- 
dim(0(p)) = pn — p{p — l)/2.) Note that the chosen quotient geometry 
differs from the one recently considered in [3], where elements of S^(p, n) are 
represented by equivalences classes ZO, O G 0(p), leading to the quotient 
representation S+(p,n) = R;''<p/0(p) . 

If ([/, R'^) e Vn,p X Pp represents A e S+(p, n), it is tempting to represent 
the tangent vectors of TaS+(p, n) by the infinitesimal variation (A, D), where 



such that Uj_ G Vn^n-p , U'^U± = 0, and Dq G Sym(p) = T/Pp. The chosen 
metric of S"'"(p, n) is merely the sum of the infinitesimal distance in Gr(p, n) 
and in Pp: 



^(t;,^2)((Ai, A), (A2, D2)) = Tr (A^^As) + A; Tr {R-^DiR-^D2R-^) , k > 0, 



generalizing (4) in a natural way. The next theorem proves that the con- 
struction endows the space S'^(p, n) with a Riemannian structure. 

Theorem 1 The space S+(p,n) = (Vn,pXPp)/0(p) endowed with the metric 
(10) is a Riemannian manifold with horizontal space 



= {(A, D) : A = U^B, B G M("-f)xf , D = RDqR, Dq G Sym(p)}. 



Furthermore, the metric is invariant with respect to orthogonal transforma- 
tions, settlings, and pseudoinversion. 

Proof: In this proof we also recap some results on quotient manifolds. We 
follow the machinery of Riemannian quotient manifold (see, e.g., [2]), except 
that we will not require that the quotient map be a Riemannian submer- 
sion. Any representative (C/, R^) lives in the structure space Vn,p x Pp. The 
tangent space to ([/, R^) is the direct sum of a vertical space and a hori- 
zontal space: T'([/^R2)(Vn,p x Pp) = V(jj^b?-^ ® 'H(u,r:^)- The equivalence class 
{(t/0, O^R^O), O G 0(p)} is called a fiber. The vertical space is the tangent 
space to the fiber at {U,R^): 



A = U^B, 
D = RDoR 



B G M^^ 



n-p)xp 



(9) 



(10) 



V(c7,fl2) = {([/n, R^VL - nR^) : n g r,o(p)}. 
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i.e., Q is a skew-symmetric matrix. The horizontal space 'H(u,r^) is by defini- 
tion complementary to V(;7r2). It is customary to represent tangent vectors 
to the quotient manifold only as elements of the horizontal space. Indeed it 
is useless to consider elements of the vertical space since they are tangent 
to the fiber, and all elements of the fiber represent the same point in the 
quotient manifold. Here the tangent space to the product manifold writes 
T(C7,H2)(V„,p X Pp) = {(S, L») : S e T,,V,,p, D e T^.Pp}. So S = A + f/Q 
with Q skew-symmetric and A = U±B. But (UQ,R^Q — QR"^) G V(^u,m)- ^'^ 
the horizontal space at T(^u^ji2^ can be chosen to be made only of the vectors 
(A, L>) given by (9): 

n^u,R^) = {(A, L>) : A = U^B, B e R(-f)xp, D = RD^R, Do G sym(p)}. 

The manifold S'^(p,n) endowed with the metric (10) is a Riemannian quo- 
tient manifold: after having chosen a Riemannian metric on the structure 
space, all wc must prove is that the induced metric on the horizontal space 
does not depend on the representative chosen. The Riemannian metric for 
the product manifold (structure space) Vn,p x Pp can be chosen as the sum of 
the natural metrics of Vn,p and Pp: S'fcfR^) 

(Xi,X2)=Tr(l]fl]2 + AfA2) + 
k Tr{DiR-'^D2R-'^), where Xi = (C/n/+ A^, A) for i = 1,2. Let Xi,X2 G 
niu,R^). We have gfJ',,,^{Xi,X2) ^ Tr {Aj A^) + k Tr {DiR'^ D^R''') . The 
first term only depends on Uj_, which is invariant along the fiber, and the 
second term is Tr (Dq, 1-^0,2)1 where Di = RDo iR, and it does not depend 
on the representative chosen. These invariances are due indeed to the invari- 
ance of the Grassman metric with respect to the representative in the Stiefel 
manifold, and the invariance properties of the natural metric on Pp. 

Concerning the last point of the theorem, an orthogonal transformation 
affects the representative (f/, i?^) by transforming U in OU with O G 0(n) 
and thus A in OA; and a scaling by transforming R^ in ij,^R^ with G 
and thus D m. ji^D. Thus they affect separately the first and second term of 
the metric, which are both invariant to each of these transformations. The 
invariance to pseudoinversion derives from the invariance to inversion of the 
natural metric on Pp, as (f/, i?~^) is a representative of a pseudoinverse of 
UR^U^ . These properties are further detailed in the proof of Theorem 2. 

The reader will note that with the proposed choice of metric and horizon- 
tal space, V([/,R2) and 'H[u,r^) are complementary, but not orthogonal. This 
is rather unusual, and it will have implications in the sequel (Proposition 1). 
However, they tend to be orthogonal as k 0. 

As a concluding remark for this section, we recall that the natural met- 
ric on the cone is not the only GL(n) -invariant metric, and a whole fam- 
ily of invariant metrics can be derived from the scalar product at identity 
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Tr(DfD2) + PTr{Di)Tr{D2), with P > (see Appendix 7.2). For in- 
stance, the metric corresponding to P = — was proposed in [20] as a 
Riemannian metric on the space of Gaussian distributions having a nonzero 
mean. This family of metrics could also easily be extended to S+(p,n). It 
should be clear that the results of Theorem 1 hold for any GL(n)-invariant 
metric on the cone and not only for the natural one. 

5 Horizontal geodesies in the structure space 
approximate geodesies in S+(p,n) 

This section provides the construction of special curves of interest connecting 

any two matrices A and B in S"'"(p, n) : horizontal geodesies in the structure 
space. First of all let us find two representatives of A and B in Vn,p x Pp 
connected by a horizontal geodesies. Let Va,Vb G Vn,p be two matrices 
that span range(A) and range(i?), respectively. The SVD of V^Va yields 
Oa, Ob e such that 

OIvJVbOb = diag((Ti, ■■■ ,ap), 1 > ai > ■ ■ ■ > > (11) 

The (Ti = cos 6i are the cosines of the principal angles < < .. < 
< 7r/2 between the two subspaces [15]. Choosing the principal vectors 
Ua = {uf, ■ ■ ■ , Up) = VaOa and Ub = {uf, ■ ■ ■ , u^) = VbOb yields a simple 
formula for the Grassman geodesic connecting range(74) and range(5) (e.g. 
[32]): 

U{t) = Ua cos (Qt) + X sin(Gt) (12) 

where O = diag(^^i, ■ ■ ■ ,6p) and X is the normalized projection of V onto the 
column space of Uj_, i.e., X = {I — UaUa)UbF where F is the pseudoinverse 
of the matrix diag(sin(6'i), ■ ■ ■ , sin(6'p)). The associated geodesic R^{t) in Pp 
must connect R\ :— U^AUa and R% :— U^BUb, that is, 

R^{t) = RAexp{t log R2^RIR^^)Ra. (13) 

Theorem 2 The singular value decomposition (11) and the geodesic curves 
(12) and (13) define a curve in S"'"(p, n) 

^A^t) = U(t)R\t)U^(t) (14) 

with the following properties: 

• 7a^s(-) connects a and B mS+(p,n), that is, ja^b{0) = A, ■jA^si^) = 
B, and -fA^Bit) G S+(p,n) Vt G [0, 1]. 
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• The curve {U (t), R'^{t)) is a horizontal lift of^A-^B{t) o-nd it is a geodesic 
in the structure space Vn^p x Pp. 

• The (squared) total length of '-yA^sif) in the Riemannian manifold 
(S+(p,n),g) is given by 

l\-iA-,B) = m\l + k \\\ogRl^RlR-^^\\l. (15) 

It is invariant with respect to pseudoinversion and to the group action 
by congruence of orthogonal transformations and scalings. 

Furthermore, the curve 7a^b(-) is uniquely defined provided that the {p — l)th 
principal angle satisfies 9p-i 7i/2. 

The proof of the theorem is given in Appendix 7.2. Viewing matrices of 
S^(p,n) as flat elhpsoids in M", the length /(7a-*b) consists of two indepen- 
dent contributions: a distance between the subspaces in which the ellipsoids 
are contained (Grassman distance), and a distance between the ellipsoids 
within a common subspace (natural distance on the cone). In this sense, 
(15) provides an exact generalization of the Riemannian distance (7) of Sec- 
tion 3. However, it is not the Riemannian distance of (S^(p,n),g) because 
ja^b{-) is not necessarily a geodesic curve (even tough it is the base curve 
of a horizontal geodesic in the structure space). 

Proposition 1 The curve 7a^b(-) is not necessarily a geodesic of{S^{p, n), g). 
Its length provides a meaningful measure of closeness between A and B which 
is not a distance because it does not satisfy the triangle inequality. 

A proof of this proposition is given by the following example: Let A — 
diag(2, 1,0) and B — diag(l,2,0). Because A and B have the same range, 
the length II'Ja^b) reduces to a distance in the cone P2 (no Grassman con- 
tribution). But A can also be connected to B via C = diag(0, 1, 2) and D = 
diag(l,0,2). The curves jA^ci'), Jc^oi-), and 7d^b(-) each involve a pure 
subspace rotation of 7r/2, which means that their total length in (S"'"(p, n), g) 
is 37r/2 (no contribution in the cone). The situation l{'yA-*B) > 37r/2 is 
clearly possible if k is large enough, showing that ja^b^-) is not necessarily 
a geodesic curve. Even if the triangular inequality is not satisfied, the pro- 
posed measure of closeness satisfies the two other properties of a distance: it 
is symmetric, and vanishes only when the two matrices coincide. 

The construction of the curve ^yA^Bi') has a geometrical meaning in 
S"'"(p,n): If a gyroscope is attached to the moving ellipsoid represented by 
U{t)R\U{t), where U{t) is given by (12), it indicates no rotation around an 
axis perpendicular to the ellipsoid during the motion < t < 1. Thus the el- 
lipsoid U{1)R\U{1) = Ub^aUb is the ellipsoid A = U{0)R\U{0) = UaRaUa 
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brought in range(5) by a rotation of minimal energy. This also justifies com- 
paring directly R\ to R% in Pp. If O 7^ / is an orthogonal matrix, it would 
violate the rotational invariance to compare R\ to O^R^^O, even though 
{UbO, O^RbO) is a valid representation of B and the distance in Grass- 
man between range(^7A) and ian.ge(UB) is unaffected by the transformation 
Ub ^ UbO. 

Even tough it is not a distance, the closeness 1{ja^b) is a meaningful 
generalization of the Riemannian distance (7) of Section 3. It reduces to 
the natural distance in the cone when A and B span the same subspace and 
it reduces to the Grassmann Riemannian distance between range(/l) and 
range (S) when A and B are rank p projectors. Furthermore, it recovers a 
maximal number of the desirable invariance properties of the natural dis- 
tance on the symmetric cone, exactly as expected from the planar example 
in Section 3. Homothetic transformations, isometrics leaving the origin fixed, 
and pseudoinversion correspond to angle-preserving transformations (when 
n > 2, these transformations are the conformal transformations of the Eu- 
clidean space) . Note we already proved in Section 3 that the invariance with 
respect to transformations which do not preserve the angles is anyway an 
impossible property to obtain, at least when p is small enough. 

The computation of the closeness (15) involves the computation of princi- 
pal angles and vectors, which is standard, and can be done via QR factoriza- 
tion at a numerical cost O(np^). The computation of the second term in (15) 
involves a symmetric generalized eigenvalue problem, which requires 0{p^) 
operations. The linear complexity in the dimension n makes the distance 
calculation efficient even in large-scale problems provided p <^n. 

The following proposition establishes a link between the closeness measure 
I and the Riemannian distance ds+(p,n) in (S'^(p, n), g^). 

Proposition 2 Let A = UaR\U^ and B = UbR%U^ be two elements of 
S"'"(p, n). We have 

< fi^A^B) - 4+(p,n)(^,^) < k max 4_^(i?|,0^i?|0). 

The right term is well defined for all O e 0{p), and it is finite since 0{p) is a 
compact set. Moreover it tends to as A; —> 0. Thus the curves defined above 
are always longer than the geodesies, but they yield good approximations of 
the geodesies when the Grassman contribution is highly penalized, i.e., k is 
small (proof in Appendix 7.2). 

We conclude this section by observing that horizontal geodesies provide 
a complete family of line search curves in S+(p,n). The curve emanating 
from A — UF^U'^ in the direction (A, D) e 'H{u,b?) admits the analytical 
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characterization 



jA{t) = U{t)R''it)U^{t) 



(16) 



where 



U{t) = UV'^ COS {rt)v + X sm{rt)V 



is the (Grassman geodesic) curve emanating from U in direction A, i.e. 
XrV = A is the compact SVD of A; and 

R\t) = Rex.p{tR~^DR-^)R 

is the (Pp geodesic) curve emanating from R^ in the direction D. An in- 
teresting feature is that neither the curve 7yi(') nor the geodesies reach the 
boundary of S^(p, n) in finite time (proof in Appendix 7.2). 

Proposition 3 For any A in S^(p, n), the curves 7a(') can be extended from 
t = — oo to t = oo. Moreover, the manifold (S"'"(p, n), gf) is geodesically 
complete. 

6 Geometric mean in S+(p,n) 

The closeness measure between A and B in the previous section provides a 
direct formula for the "halfway" matrix 



where K = Ra{R']^R\R']^Y^'^Ra is the Riemannian mean of i?^ and R\ in 
Pp, W = cos(r/2)f/A + sin(r/2)X is the Riemannian mean of range(74) and 
range(5), and where Ra, Rb,Ua, X are defined via the SVD (11). 

The mean is uniquely defined in the generic case 9p ^ as the mean 
in Grassman between range(A) and range(S) is uniquely defined. When 
6p = n/2 and 9p-i ^ t^I'^-, according to Theorem 2 1{'^a-^b) is uniquely 
defined, but there are two Riemannian means of range(A) and range(i?) in 
Grassman. So the mean (17) is well-defined but there are two midpoints. If 
9p-i = 7r/2, there is an infinity of midpoints. Indeed, let r be the number of 
principal angles equal to 7r/2. Looking at the proof of Theorem 2 in Appendix 
7.2, we see that the midpoint will necessarily be a matrix of the form (17) 
where Ua, Ra are replaced by f/^P, P^RaP, with P G 0(p) representing an 
arbitrary rotation on the span of the p — r last principal vectors. From now 
on we will systematically assume that 6p ^ n/2. 



AoB^ WKW^, 



(17) 
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The geometric mean A o B is closely related to the popular geometric 
mean A^B of Ando [5, 29, 4]. One definition of Ando mean is to consider 
A^B as the solution of the extremal problem 

max {xhO : ol . (18) 



BJ - 

For full-rank positive definite matrices A and B, the solution of (18) is given 
by X = (^-1/2^^-1/2^ 1/2^1/2^ ^Yi&t is, Ando mean coincides with the 
geometric mean (17). In contrast, for matrices A and B in S"'"(p, n), the 
geometric mean (17) differs from the Ando mean (18). In particular the 
geometric mean is by definition rank-preserving whereas the solution of (18) 
has a rank which is upper bounded by dim(range(74) fl range(i?)). A simple 
example is provided by the example already discussed in Section 3. The 

geometric mean of /l=diag(4,0) and i?=diag(0,l) is A o B = ^ ^ while 

the solution of (18) is A^B — • This is easy to check because the 

Ando mean can be obtained by density: letting A' — A + el and B' — B + el 
the mean corresponds to the limit of A'V2(^/-i/2^/^/-i/2^ 1/2^/1/2 ^ ^ 0. 

Apart from this important difference, both the geometric mean (17) and 
the Ando mean enjoy most desirable properties of a matrix geometric mean 
[4], listed in the proposition below. 

Proposition 4 The mean Ao B possesses the properties listed below. 

i. Joint homogeneity aA o pB = (a^Y^^A o B. 
a. Permutation invariance A o B — B o A. 

Hi. Monotonicity. li A < Aq (i.e. (Aq — A) is a positive matrix) and 
B < Bq, the means are comparable and verify Ao B < Aqo Bq. 

iv. Congruence invariance. For any (/x, P) e ]RxO(n) we have {jiP'^ AiJiP)o 
(plP^BhP) = nP^{A o B)ijlP. 

V. Self-duality {A o By = {A^ o B"^), where "f" denotes pseudoinversion. 

The mean (17) may prove useful to generalize to low-rank positive semidef- 
inite matrices the growing use of the geometric mean in applications re- 
quiring interpolation and the filtering of positive definite matrices; see e.g. 
[26, 6, 7, 34]. 
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A particular case of interest is the set of rank p projectors 

{P e M"^7 = P, P^ = P, Tr (P) = p}, (19) 

which is in bijection with the Grassman manifold of p-dimensional subspaces. 
Not surprisingly, the geometric mean (17) of A and B in (19) agrees with 
the Riemannian mean of range(yl) and rangc(i?) on Gr(j9, n). On this set, 
the geometric mean also agrees with the Riemannian mean in S^(p, n) since 
ds+{p,n){A,B) = (iGr(p,n)(range(y4),range(5)) for A, B projectors (see Ap- 
pendix 7.2). 

On the other hand, when A and B have the same range, the geomet- 
ric mean (17) docs not necessarily agree with the Riemannian mean in 
S"'"(p, n) since the connecting curve 7a^b(-) is not necessarily a geodesic. 
Note however that the geometric mean is a good approximation of the Rie- 
mannian mean when k is small. In fact the geometric mean (17) is a 
more natural extension of the Ando mean than the Riemannian mean since 
AoB^ U{R\#Rl)U^ for A = UR\U^ and B = URlU^. 

7 Conclusion 

This paper generalizes the Riemannian geometry of the symmetric cone Pn to 
the manifold S"'"(p,n) of positive semidefinite matrices of fixed-rank p. The 
generalization is based on the quotient geometry S+(p, n) = (Vn,p x Pp) /0(p) 
that leads to a natural metric with decoupled contribution in Gr(p, n) and the 
cone Pp. This geometry leads to an explicit and natural notion of closeness 
in S"'"(p,n), which in turn provides a computable and natural definition of 
Riemannian distance and rank-preserving geometric mean. 

The proposed computational tools may prove useful in applications in- 
volving computations with low-rank approximations of large-scale positive 
definite matrices. Such applications have already appeared in MRI tensor 
computing [26, 25, 14, 6], and in radar processing [7]. Particular areas where 
the tools may prove useful in the future are the growing use of kernel-based 
methods, and low-rank approximations in machine learning and in bioinfor- 
matics. Kernel learning [33] , kernel completion [30] , and the use of Bregman 
divergence to address matrix nearness problems [11, 19] (the natural distance 
on the cone is a Bregman divergence) are exemplary illustrations of areas that 
could benefit from the computational tools introduced in this paper. 
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Appendix 

7.1 Reductive homogeneous space and invariant met- 
rics 

Let us choose the identification P„ = Gl(n)/0(n) as a starting point. Ev- 
ery matrix of the Lie algebra gl{n) = T/Gl(n) = R"^" can be expressed as 
the sum of its symmetric part and its skew-symmetric part. Thus Qi{n) = 
m -|- so(n) is the direct sum of the symmetric matrices m and the Lie alge- 
bra so(n) = T/0(n), made of skew-symmetric matrices. Moreover, for any 
5 e m and O e 0(n) we have Ado{S) = OSO'^ = OSO^ e m. Thus 
^c^o(n)(Tn) C m. The existence of m such that Ql{n) = m + so{n) is a direct 
sum, and A(io(n)(nT) C m proves that Gl(n)/0(n) can be called a reductive 
homogeneous space. 

We provide some more information for the readers who are familiar with 
standard ideas of Lie group theory. Let {Z{t), t > 0} be a trajectory in 
Gl(n). {Zi{t) = Z{t)0{t), t > 0} with 0{t) e 0(n) represents the same 
trajectory in P„. Let uj{t) = Z~'^^Z G gl{n) and uji{t) = Z^^^Zi. We have 
u!{t) = Ado{t){i-Oi{t) — 0{t)~^-^0{t)) (see [18]). Decomposing u!,u!i in their 
symmetric and skew-symmetric parts, we see that any A(io(n)-invariant scalar 
product on m will allow the construction of a well-defined GL(n)-invariant 
metric on the homogeneous space Gl(n)/0(n) corresponding to this scalar 
product at /. All /l(io(n)-invariant scalar products on the symmetric matrices 
are given up to a constant scale factor by Tr (^XfX2) + /3Tr (Xi) Tr {X2) with 
f3 > —K Indeed, they are derived from rotationally invariant norms 
on symmetric matrices. Thus they can only depend on the scalar invariants 
Tr(X), Tr(X^), Tr(X^), etc. As they are quadratic functions they can 
only depend on Tr (X)^ and Tr(X^). The condition on /? ensures positive 
definiteness (see, e.g., [27]). 

7.2 Several proofs 



Proof of Theorem 2 Let U{t) be a geodesic linking Ua and in Vn,p 
and R'^{t) a geodesic hnking R\ and R% in Pp, where {Ua, Ra) and {Ub, R%) 
are representatives of A and B defined via SVD (11). By corollary 3.57 of 

[24] {U (t), R'^{t)) is a geodesic curve in the structure space Vn,p x Pp endowed 
with the metric g^^. The choice of the representatives is such that for every 
t > 0, -^U (t) is orthogonal to U (t), as it is an element of the horizontal space 
of Gr(p, n) = Vn,p/0(p). This has two consequences. First the length of the 
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geodesic in the structure space is 

(f{A, B) = 4,(p^„)(range(^), range(5)) + k 4^{R\, RI). 

Because the Riemannian distance in Gr(p, n) is ||©||f [12, 1] and the Rie- 
mannian distance in Pp is ||logi?^^i?^i?^^||i^, we obtain 

d\A,B) = ||e||^ + fc WlogR^'RlR^'Wl. 

Then ^{U{t), R^{t)) belongs to the horizontal space 'H(u{t),B?{t)) for all t > 0, 
that is, {U{t),R^{t)) is a horizontal curve. As g^^ induces the metric (10) 
on the quotient manifold, the length I^ja^b) of Ja^b in the quotient space 
(Vn,p X Pp)/0(p) is also d{A,B), proving (15). 

Uniqueness. First of all suppose the crj's are distinct. According to the 
uniqueness of the SVD in (11), the matrices Ua and Ub are unique (up to 
a joint multiplication of any columns by —1, which is an orthogonal trans- 
formation) and do not depend on the choice of Va and Vb- Under the less 
restrictive assumption crp_i > 0, suppose there exists 1 < ^ < p — 1 such 
that CTj = (Tj+i. Then cTj > and the SVD yields uf,uf,u^'j^,uf_^i such 
that {ufYuf = {u^i)'^uf_^i = (Ti. For the sake of simplicity we assume 
there are only 2 principal vectors associated to (jj. The generalization to 
an arbitrary number of vectors is straightforward and leads to the same 
conclusion. The SVD yields non-unique principal vectors, since any nor- 
malized linear combination of Ui and i^i+i is still a principal vector. Let 

= auf + hufj^i and = cuf + duf^^^ be other principal vectors asso- 
ciated to (Tj. By definition {u'^)^u^ = ai and {u'^Y'u'^ = {u^)^u^ = 1. 
Since (-uf )"^-uf^]^ = we have (ac + bd)ai = ai. It implies a = c and 
b — d (Cauchy-Schwartz equality) and, necessarily, and are obtained 
from uf,u^i, and uf,uf^i via the same orthogonal transformation (unless 
(Tj = 0, which is impossible). We proved the principal vectors Ua and Ub 
are defined up to a joint multiplication on the right by a block orthogonal 
matrix P G 0(p), each block being associated to the same eigenvalue. Let 
UaP = U'a and UbP = U'^, and let A = U'aR'HU'a^ and B = U'bR'^{Ub)'^ ■ 
We have R'l = P^RaP, R'i = P^R^P and /(ta-^b) is unchanged since 
dpp(it!^, it!|) = dp^{R'A, R'q). Note that the reciprocal of this result is also 
true. 

Invariances. Using the preceding paragraph, we know that (15) is well 
defined under the basic assumption crp_i > 44> 9p-i < 7r/2. Thus the 
following features are sufficient to complete the proof: 

• For /X G M we have jj.'^A = Ua{i^'^R^a)Ua- Since we have dp^{ii^R\, n'^R%) = 
dp^{R%Rl) we have lij^2A^^2B) = 1{ia^b)- 



METRICS FOR POSITIVE SEMIDEFINITE MATRICES 



20 



• Let O e 0(n), and let Ua and Ub he the principal vectors associ- 
ated to A, B. Then f/^ = OUa and U'^ = OUb are principal vectors 
associated to OAO^ and OBO^ since orthogonal transformations pre- 
serve the angles: {OUaY^OUb = U]^Ub, so the Grassman distance 
is unchanged. Since OAO^ = U'^R\U'^ and OBO^ = U'^RlU'^ 
we see that R\, R% are also unchanged by the transformation and 
Kioao^-*obot) = /(ta-^b)- 

• A^ — UaRa^U^ is the pseudoinverse of A. 1{^a^b) — /(7At^_Bt) = 
^Pp(-^A > -^B ) - dp^(R% R%) = 0. 

Proof of Proposition 2 First, as 7a^b induces a non-minimal path con- 
necting A and B in S+(p,n) we have /^(7a-+b) > ds+{p,n){A, B). Then, let 
7(t) be a geodesic linking A and i? in S"'"(p, n), with 7(0) = A, 7(1) = B 
and (is+(p,n)(^, -B) its length. As proved in the sequel 7(t) can be glob- 
ally horizontally hfted in the structure space Vn,p x Pp. The horizontal lift 
7(t) at {Ua, Ra) is a horizontal curve in the structure space Vn,p x Pp con- 
necting the two fibers. Therefore there exists O G 0{p) such that 7(1) = 
{UbO, O^RlO), where {Ub, Rl) is defined via the SVD (11). The length of 7 
in the structure space is also (is+(p,n)(^) B). Since 7 is not necessarily a min- 
imal curve in the structure space we have 'i|+(p n) B) > dy^^{UA, UbO) + 
k dp {R\,O^R%0). Because distances in Gr(p, n) are shorter than those in 
Vn,p,'this implies 4+(p,n)(^, B) > 4^(p^)(range(C/^), range(C/BO))+A: dl^{R\, C^RlO). 
Finally we proved 

1\ia^b) < 4r(p,n)(range([/A), range([/B)) + k dl^{R\, 0^i?|0) + k 4^(0^i?|0, 

< dl,^^^^){A,B) + kdl^{Rl,O^RlO) 

Proof of Proposition 3 The structure space is geodesically complete 
(it is a product of two geodesically complete Riemannian manifolds whose 
geodesies are product of geodesies). Thus horizontal geodesies are defined 
for all t e M. 

A manifold is geodesically complete if and only if it is complete as a 
metric space (see Theorem 10.3 Ch I in [16]). Thus the structure space is 
complete. To prove that the quotient space S^(p, n) is geodesically complete 
it suffices to prove that it is complete. Let {um} be a Cauchy sequence in the 
quotient space S+(p,n). One can find a broken geodesic 7(t) linking the Wj's 
such that 7(ti) = Ui. Consider an horizontal hft 7 of 7 (it is well defined as 
proved below). The distance between ^{U) and ^{tj) in the structure space is 
ds+{v,r,){l{U),l{tj)) = ds+(p,n)(iij,iij) and thus {^{tm)} is a cauchy sequence 
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in the structure space. Thus it converges to a point p — {U,R^). Let us 
prove {ujn\ converges io p = UR^U^ . Thanks to the exponential map in the 
structure space, 7(ti) can be identified to an element = (Aj + UQi,Di) G 
^(c/,i?2)Vn,p X Pp for i large enough. The geodesic distance in the structure 
space between ^{U) and p is 5-^^(6,6) = W^if + WW + ^ Tr ((A^^^)^) 
and thus \\^i\\ 0. Thus the coordinate of ^{U) along the fiber tends to 
zero so ds+ (p,n) {Ui Rj , UR^U^) 0. 

Any geodesic 7 of S^(p,n) can be globally horizontally lifted. This re- 
sult was shown by R. Hermann when the horizontal and vertical spaces are 
orthogonal. The result holds in our case and the proof is unchanged. Let us 
recall the main points of the proof : let 7^ is the maximal horizontal lift of 
7 at p, defined on, say [0, e[. It is enough to show that 7p can be extended 
at e. The continuity of 7 and the completeness of the structure space imply 
the existence of a limit point q = limt^^ 7p(^)- 
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